2-2

CompConj[A_] := A /.Complex[0,n_]->-Complex[0,n]

Quantum Vortices in a Square Well

Nicholas Wheeler
Reed College Physics
2 May 2008

Introduction

I was asked recently by the acting editor of AJP to referee a paper "Quantum-mechanical phase vortices in the two-dimensional infinite square well" by Meredith A. Mahoney, David M. Paganin & Helen M. L. Faulkner. I made heavy use of Mathematica in preparing and illustrating my report (which survives as a couple of pdf files), but the notebooks generated by that effort were considered to be mere worksheets, and were left in a very disorganized state. I was motivated by MPF to undertake the work recorded in "Quantum vortices in a circular well," which I had imagined would be relatively more interesting/instructive, but were discovered to be less so. So I have come belatedly to value the MPF work for its computational tractability and the vividness of the results to which it leads, and am motivated to undertake this more careful account of my own work in this subject area.

Some preliminary remarks of a general nature

Let  Ψ[x, y] be written in polar form

Ψ=R Exp[(I/ℏ)S];
Ψstar=R Exp[-(I/ℏ)S];

with

R=r[x,y,t];
S=s[x,y,t];

Supposing Ψ[x, y] to be a solution of the free-particle Schrödinger equation, we construct

Schrödinger=Exp[-(I/ℏ)S](-((ℏ^2)/(2m))(D[Ψ,{x,2}]+D[Ψ,{y,2}])-I ℏ D[Ψ,t])//Expand

from which we extract and massage the real and imaginary parts (which must vanish separately):

((CompConj[Schrödinger]+Schrödinger)/(2  r[x,y,t]))//Expand
((CompConj[Schrödinger]-Schrödinger)/(2 ℏ I))//Expand

The first equation can be written

Derivative[0,0,1][s][x,y,t]+(((Derivative[0,1,0][s][x,y,t])^2)/(2 m))+(((Derivative[1,0,0][s][x,y,t])^2)/(2 m))=(((ℏ^2)( Derivative[2,0,0][r][x,y,t]+Derivative[0,2,0][r][x,y,t]))/(2 m r[x,y,t]))

where in classical theory the expression on the left would vanish by appeal to Hamilton-Jacobi theory. The second equation finds employment at the end of the next paragraph:

The probability density is given then by

Ψstar Ψ

and the probability current by

(Subscript[J, x])=(ℏ/(2m I))(Ψstar D[Ψ,x]-Ψ D[Ψstar,x])//Simplify
(Subscript[J, y])=(ℏ/(2m I))(Ψstar D[Ψ,y]-Ψ D[Ψstar,y])//Simplify

The continuity equation that provides a local account of probability conservation lends interest to a construction

D[((r[x,y,t])^2),t]+D[((((r[x,y,t])^2) Derivative[1,0,0][s][x,y,t])/m),x]+D[((((r[x,y,t])^2) Derivative[0,1,0][s][x,y,t])/m),y]//Expand

that after only a little massaging leads to a construction that has already been shown to vanish;

((%)/(r[x,y,t]))//Expand

%==((CompConj[Schrödinger]-Schrödinger)/(2 ℏ I))//Expand

The vector-valued construction

((probability current)/(probability density))

has the physical dimension of a "velocity." It was shown above that

quantum mechanical velocity field ≡ ((probability currrent)/(probability density))=(1/m)Grad[S]

Exists except

The velocity field

(({{(Subscript[V, x])[x,y,t]}, {(Subscript[V, y])[x,y,t]}}))=(((({{(Subscript[J, x])[x,y,t]}, {(Subscript[J, y])[x,y,t]}})))/(P[x,y,t]))=(1/m)Grad[S]

exists except where the probability density P[x,y,z,t] vanishes (which is to say: where Ψ[x,y,z,t] vanishes), and might exist even at such points if  J[x,y,z,t] happened also to vanish in just the right way. At those points the phase must be undifferentiable (discontinuous).

In fluid dynamics (see §1.2—"Rotation and Vorticity"—pages 19-31 in A. J. Chorin & J. E. Marsden, A Mathematical Introduction to Fluid Mechanics (2nd edition, 1990)) the vector-valued field produced by constructing the curl of a velocity field is called the "vorticity":

Ω=Curl[V]

For the "quantum mechanical probability fluid" one has

Ω=(1/m)Curl[ Grad[S]]=0

except at singular points (points where phase discontinuities cause  Grad[S]  to be undefined).

Quantum vorticity of a toy model

Let us agree, now and hencforth, to adopt units in which

2m = ℏ = 1

The 2-dimensional free-particle Schrödinger equation then reads

-(((∂_x)^2)+((∂_y)^2))ψ=i (∂_t)ψ

I have need now of this

A seldom remarked point: Let

Ψ[x,y,t,α,β]

be a solution of the Schrödinger equation that depends upon certain parameters {α, β, ...}. Then (by linearity)

(Subscript[Ψ, mn])[x,y,t,α,β]=((((∂_α)))^m)((((∂_β)))^n)Ψ[x,y,t,α,β]

is also a solution. So are all real/complex linear combinations of such functions. Note, however, that functions produced in this way will, in general, not mimic the boundary conditions satisfied by the mother function Ψ. Nor will they, in general, be normalized.

Look, for example, to the function

Ψ=Exp[I(α x+β y-((α^2)+(β^2)) t)]

which clearly does satisfy the free-particle Schrödinger equation:

-D[Ψ,{x,2}]-D[Ψ,{y,2}]-I D[Ψ,{t,1}]//Simplify

Construct

(Subscript[Ψ, 11])=-I D[Ψ,{α,1}]+ D[Ψ,{β,1}]//Simplify

which we (by the "seldom remarked point") know, but will anyway demonstrate, satisfies the Schrödinger equation:

-D[(Subscript[Ψ, 11]),{x,2}]-D[(Subscript[Ψ, 11]),{y,2}]-I D[(Subscript[Ψ, 11]),{t,1}]//Simplify

From

(Subscript[Ψ, 11])CompConj[(Subscript[Ψ, 11])]//ComplexExpand

we are led to define the density

P[x_,y_,t_,α_,β_]:=(x^2)+(y^2)-4 t x α+4 (t^2) (α^2)-4 t y β+4 (t^2) (β^2)

which, though everywhere non-negative, is not normalized, so cannot be interpreted as "probability density." Introducing the associated current (which is not "probability current")

-I(CompConj[(Subscript[Ψ, 11])] D[(Subscript[Ψ, 11]),x]-(Subscript[Ψ, 11]) D[ CompConj[(Subscript[Ψ, 11])],x])//Expand
-I(CompConj[(Subscript[Ψ, 11])] D[(Subscript[Ψ, 11]),y]-(Subscript[Ψ, 11]) D[ CompConj[(Subscript[Ψ, 11])],y])//Expand

JX[x_,y_,t_,α_,β_]:=-2 y+2 (x^2) α+2 (y^2) α-8 t x (α^2)+8 (t^2) (α^3)+4 t β-8 t y α β+8 (t^2) α (β^2)
JY[x_,y_,t_,α_,β_]:=2 x-4 t α+2 (x^2) β+2 (y^2) β-8 t x α β+8 (t^2) (α^2) β-8 t y (β^2)+8 (t^2) (β^3)

we confirm that we do have the anticipated local conservation law (continuity equation)

D[P[x,y,t,α,β],t]+D[JX[x,y,t,α,β],x]+D[JY[x,y,t,α,β],y]//Simplify

Needs["VectorFieldPlots`"]

Manipulate[Plot3D[P[x,y,t,α,β],{x,-1,1},{y,-1,1}, PlotRange->{0,20}, Ticks->False, PlotLabel->"Probability Density", LabelStyle->Directive[Bold, Red]],{t,0,1},{α,-1,1},{β,-1,1}]
Manipulate[VectorFieldPlot[{JX[x,y,t,α,β],JY[x,y,t,α,β]},{x,-1,1},{y,-1,1}, PlotLabel->"Probability Current", LabelStyle->Directive[Bold, Red]],{t,0,1},{α,-1,1},{β,-1,1}]

The figure reveals a moving swirl when α and β are assigned non-zero values.

The probability density vanishes if and only the wave function itself vanishes, which requires that its real and imaginary parts both vanish. From

(Subscript[Ψ, 11])

we learn that the probability density vanishes (and the velocity field becomes singular) at a point  {x,y} that moves linearly

x-2α t=0
y-2β t=0

P[x,y,t,α,β]/.{x->2α t, y->2β t}

The velocity field is defined

VX[x_,y_,t_ ,α_,β_]:=((JX[x,y,t,α,β])/(P[x,y,t,α,β]))
VY[x_,y_,t_ ,α_,β_]:=((JY[x,y,t,α,β])/(P[x,y,t,α,β]))

VX[x,y,t,α,β]

but to humor Mathematica—to keep her from complaining

VX[x,y,t,α,β]/.{x->2α t, y->2β t}

—we adopt this modified definition:

Clear[α]

VVX[x,y,t,α,β]

VVX[x_,y_,t_ ,α_,β_]:=((JX[x,y,t,α,β])/(P[x,y,t,α,β]+0.00000001))
VVY[x_,y_,t_ ,α_,β_]:=((JY[x,y,t,α,β])/(P[x,y,t,α,β]+0.00000001))

Manipulate[VectorFieldPlot[{VVX[x,y,t,α,β],VVY[x,y,t,α,β]},{x,-1,1},{y,-1,1}, PlotLabel->"Velocity Field", LabelStyle->Directive[Bold, Red]],{t,0,1},{α,-1,1},{β,-1,1}]

The vorticity (curl of the velocity field) is in 2-dimensional theory a pseudoscalar field, known from preceding remarks to vanish

D[VY[x,y,t,α,β],x]- D[VX[x,y,t,α,β],y]//Simplify

except at the singular point. To expose the singularity to grapic scrutiny, we proceed from

D[VVY[x,y,t,α,β],x]- D[VVX[x,y,t,α,β],y]//Simplify

to the definition

ΩΩ[x_,y_,t_ ,α_,β_]:=((4.`*^-8-4.`*^-8 y α+4.`*^-8 x β)/(((1.`*^-8+(x^2)+(y^2)-4.` t x α+4.` (t^2) (α^2)-4.` t y β+4.` (t^2) (β^2)))^2))

This function has a large value at the singularity (located at the origin of the {x,y}-plane when t = α = β =0)

ΩΩ[0,0,0,0,0]

but falls to very small values at close-by points

ΩΩ[0.01,0.01,0,0,0]
ΩΩ[0.1,0.1,0,0,0]

As previously remarked, we expect the singular point to be a point at which the phase of the wave function

ΨΨ[x_,y_,t_,α_,β_]:=(E^(I (x α+y β-t ((α^2)+(β^2))))) (x+I (y+2 I t α-2 t β))

becomes undifferentiable. The close relationship between the vorticity and phase singularities is illustrated in the following sequence of paired figures

Plot3D[ΩΩ[x,y,0,0,0],{x,-1,1},{y,-1,1}, PlotRange->{-100,100}, Ticks->False, PlotPoints->100, PlotLabel->"Vorticity", LabelStyle->Directive[Bold, Red]]
Plot3D[
Arg[ΨΨ[x,y,0,0,0]],{x,-1,1},{y,-1,1}, PlotRange->{-π,π}, Ticks->{{},{},{-π,π}}, PlotPoints->100,PlotLabel->Style["Phase Singularity",Bold,Red]]

Plot3D[ΩΩ[x,y,.25,.25,.25],{x,-1,1},{y,-1,1}, PlotRange->{-100,100}, Ticks->False, PlotPoints->100, PlotLabel->"Vorticity", LabelStyle->Directive[Bold, Red]]
Plot3D[
Arg[ΨΨ[x,y,.25,.25,.25]],{x,-1,1},{y,-1,1}, PlotRange->{-π,π}, PlotPoints->100, Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase Singularity",Bold,Red]]

Plot3D[ΩΩ[x,y,.5,.5,.5],{x,-1,1},{y,-1,1}, PlotRange->{-100,100}, Ticks->False, PlotPoints->100, PlotLabel->"Vorticity", LabelStyle->Directive[Bold, Red]]
Plot3D[
Arg[ΨΨ[x,y,.5,.5,.5]],{x,-1,1},{y,-1,1}, PlotRange->{-π,π},  PlotPoints->100, Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase Singularity",Bold,Red]]

To expose most clearly the anatomy of the phase singularity one can (as MPF advocate) construct the first-order Taylor expansion of  ΨΨ[x,y,t,α,β] about the singular point x = 2 α t, y = 2 β t.

Series[ΨΨ[2α t+dx,2β t+dy,t,α,β],{dx,0,1},{dy,0,1}]

Use that information to construct the following function, which is "tangent to the wave function at the point where it vanishes"

Clear[ψ]

ψ[dx_,dy_,t_,α_,β_]:=(E^(I t (α^2)+I t (β^2)))(dx+I dy)

and then examine the phase of that relatively simple function:

Plot3D[
Arg[ψ[dx,dy,.5,.5,.5]],{dx,-.1,.1},{dy,-.1,.1}, PlotRange->{-π,π}, PlotPoints->100, Ticks->{{},{},{-π,π}},PlotLabel->Style["Close-up  Detail of the Phase Singularity",Bold,Red]]

Energy eigenfunctions for the one-dimensional square well of infinite depth

We seek normalized solutions of the time-independent Schrödinger equation

-((∂_x)^2)ψ[x]=ℰ  ψ[x]

that conform to the boundary conditions

ψ[0]=ψ[1]=0

We are led to the functions

f1[x_]:=(Sqrt[2])Sin[1π x]
f2[x_]:=(Sqrt[2])Sin[2π x]
f3[x_]:=(Sqrt[2])Sin[3π x]

and their higher order relatives.

(∫_0^1)((f1[x])^2)\[DifferentialD]x
(∫_0^1)((f2[x])^2)\[DifferentialD]x
(∫_0^1)((f3[x])^2)\[DifferentialD]x

((-D[f1[x],{x,2}])/(f1[x]))
((-D[f2[x],{x,2}])/(f2[x]))
((-D[f3[x],{x,2}])/(f3[x]))

Generally, (Subscript[ℰ, n])=(n^2)(π^2) so the buzzing eigenfunctions (solutions of the time-dependent Schrödinger equation) are

F1[x_,t_]:=f1[x]Exp[-I (1^2)(π^2)t]
F2[x_,t_]:=f2[x]Exp[-I (2^2)(π^2)t]
F3[x_,t_]:=f3[x]Exp[-I (3^2)(π^2)t]

-D[F1[x,t],{x,2}]-I D[F1[x,t],t]
-D[F2[x,t],{x,2}]-I D[F2[x,t],t]
-D[F3[x,t],{x,2}]-I D[F3[x,t],t]

Energy eigenfunctions for the two-dimensional square well of infinite depth

The system separates into two copies of the system just studied, so we are led immediately to the buzzing normalized eigenfunctions

F11[x_,y_,t_]:=f1[x]f1[y]Exp[-I (π^2)((1^2)+(1^2))t]
F12[x_,y_,t_]:=f1[x]f2[y]Exp[-I (π^2)((1^2)+(2^2))t]
F13[x_,y_,t_]:=f1[x]f3[y]Exp[-I (π^2)((1^2)+(3^2))t]
F21[x_,y_,t_]:=f2[x]f1[y]Exp[-I (π^2)((2^2)+(1^2))t]
F22[x_,y_,t_]:=f2[x]f2[y]Exp[-I (π^2)((2^2)+(2^2))t]
F23[x_,y_,t_]:=f2[x]f3[y]Exp[-I (π^2)((2^2)+(3^2))t]
F31[x_,y_,t_]:=f3[x]f1[y]Exp[-I (π^2)((3^2)+(1^2))t]
F32[x_,y_,t_]:=f3[x]f2[y]Exp[-I (π^2)((3^2)+(2^2))t]
F33[x_,y_,t_]:=f3[x]f3[y]Exp[-I (π^2)((3^2)+(3^2))t]

These functions are obviously/demonstrably normalized

(∫_0^1)((∫_0^1)F12[x,y,t]CompConj[F12[x,y,t]]\[DifferentialD]x)\[DifferentialD]y

and obviously/demonstrably satisfy the Schrödinger equation:

-D[F12[x,y,t],{x,2}]-D[F12[x,y,t],{y,2}]==I D[F12[x,y,t],t]

Their probability densities are (of course) time-independent

Prob1212[x_,y_,t_]:=F12[x,y,t]CompConj[F12[x,y,t]]

Prob1212[x,y,t]

Plot3D[4 ((Sin[π x])^2) ((Sin[2 π y])^2),{x,0,1},{y,0,1},Ticks->False,PlotLabel->"Static Probability Density of a Buzzing Eigenfunction", LabelStyle->Directive[Bold, Red]]

and the associated probability currents vanish:

-I(CompConj[F12[x,y,t]]D[F12[x,y,t],x]-F12[x,y,t]D[CompConj[F12[x,y,t]],x])
-I(CompConj[F12[x,y,t]]D[F12[x,y,t],y]-F12[x,y,t]D[CompConj[F12[x,y,t]],y])

To produce motion—which in quantum mechanics is an interference effect—one must take such functions in (normalized) linear combination. In the simplest instance we have constructs of the form

Clear[Ψ,μ,ϕ]

Ψ[x_,y_,t_,μ_,ϕ_]:=Cos[μ]F21[x,y,t]+(E^(I ϕ))Sin[μ]F12[x,y,t]

Ψ[x,y,t,μ,ϕ]

where μ is a "mixing angle" and ϕ a "relative phase" factor. For no better reason than to achieve agreement with MPF I set

FIRST EXAMPLE

For no better reason than to achieve agreement with MPF, I set

μ=ArcCos[(2/3)]
ϕ=ArcTan[2]

Cos[μ]
(E^(I ϕ))Sin[μ]//ComplexExpand

Ψ=(2/3)F21[x,y,t]+((1+2I)/3)F12[x,y,t]

ψ=Ψ/.t->0

ψstar=(4/3) Sin[2 π x] Sin[π y]+((2/3)-((4 I)/3)) Sin[π x] Sin[2 π y]

ψstar ψ//ComplexExpand

P[x_,y_]:=(16/9) ((Sin[2 π x])^2) ((Sin[π y])^2)+(16/9) Sin[π x] Sin[2 π x] Sin[π y] Sin[2 π y]+(20/9) ((Sin[π x])^2) ((Sin[2 π y])^2)

(∫_0^1)((∫_0^1)P[x,y]\[DifferentialD]x)\[DifferentialD]y

Plot3D[P[x,y],{x,0,1},{y,0,1},Ticks->False,PlotLabel->"Initial Probability Density: Example 1", LabelStyle->Directive[Bold, Red]]

It is evident from the figure, but I now verify, that the initial probability density vanishes at the middle of the box:

P[1/2,1/2]

-I(ψstar D[ψ,x]-ψ D[ψstar,x])//ComplexExpand
-I(ψstar D[ψ,y]-ψ D[ψstar,y])//ComplexExpand

JX[x_,y_]:=-(64/9) π Cos[2 π x] Sin[π x] Sin[π y] Sin[2 π y]+(32/9) π Cos[π x] Sin[2 π x] Sin[π y] Sin[2 π y]
JY[x_,y_]:=(64/9) π Cos[2 π y] Sin[π x] Sin[2 π x] Sin[π y]-(32/9) π Cos[π y] Sin[π x] Sin[2 π x] Sin[2 π y]

VectorFieldPlot[{JX[x,y],JY[x,y]},{x,0,1},{y,0,1},PlotLabel->"Initial Probability Current: Example 1", LabelStyle->Directive[Bold, Red]]

((JX[x,y])/(P[x,y]+0.000001))

VVX[x_,y_]:=((JX[x,y])/(P[x,y]+0.000001))
VVY[x_,y_]:=((JY[x,y])/(P[x,y]+0.000001))

VectorFieldPlot[{VVX[x,y],VVY[x,y]},{x,0,1},{y,0,1},PlotLabel->"Initial Velocity Field: Example 1", LabelStyle->Directive[Bold, Red]]

D[VVY[x,y],x]-D[VVX[x,y],y]//Simplify

Ω[x_,y_]:=(0.00028073541407543063` Cos[2 π x] Cos[2 π y] Sin[π x] Sin[π y]-0.00007018385351885766` Cos[π x] Cos[π y] Sin[2 π x] Sin[2 π y])/(((1.`*^-6+1.7777777777777777` ((Sin[2 π x])^2) ((Sin[π y])^2)+1.7777777777777777` Sin[π x] Sin[2 π x] Sin[π y] Sin[2 π y]+2.2222222222222223` ((Sin[π x])^2) ((Sin[2 π y])^2)))^2)

Plot3D[Ω[x,y],{x,0,1},{y,0,1}, PlotRange->{0,10},PlotPoints->{100}, Ticks->False,PlotLabel->"Initial Vorticity: Example 1", LabelStyle->Directive[Bold, Red]]

Plot3D[Arg[(4/3) Sin[2 π x] Sin[π y]+((2/3)+((4 I)/3)) Sin[π x] Sin[2 π y]],{x,0,1},{y,0,1}, PlotRange->{-π,π},PlotPoints->{100},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Initial Phase: Example 1",Bold,Red]]

SECOND EXAMPLE

Here, following again in the steps of MPF, I use

Cos[(π/4)]
(E^(I π/2))Sin[(π/4)]

to assemble

ψ2[x_,y_,t_]:=(1/(Sqrt[2]))F31[x,y,t]+(I/(Sqrt[2]))F13[x,y,t]

ψ2[x,y,t]

ψ2star[x_,y_,t_]:=(Sqrt[2]) (E^(+10 I (π^2) t)) Sin[3 π x] Sin[π y]-I (Sqrt[2]) (E^(+10 I (π^2) t)) Sin[π x] Sin[3 π y]

Proceeding now through the familiar sequence of constructions…

ψ2star[x,y,0]ψ2[x,y,0]//ComplexExpand

P2[x_,y_]:=2 ((Sin[3 π x])^2) ((Sin[π y])^2)+2 ((Sin[π x])^2) ((Sin[3 π y])^2)

Plot3D[P2[x,y],{x,0,1},{y,0,1}, Ticks->False,PlotLabel->"Initial Probability Density: Example 2", LabelStyle->Directive[Bold, Red]]

-I(ψ2star[x,y,0] D[ψ2[x,y,0],x]-ψ2[x,y,0] D[ψ2star[x,y,0],x])//ComplexExpand
-I(ψ2star[x,y,0] D[ψ2[x,y,0],y]-ψ2[x,y,0] D[ψ2star[x,y,0],y])//ComplexExpand

JX2[x_,y_]:=-12 π Cos[3 π x] Sin[π x] Sin[π y] Sin[3 π y]+4 π Cos[π x] Sin[3 π x] Sin[π y] Sin[3 π y]
JY2[x_,y_]:=12 π Cos[3 π y] Sin[π x] Sin[3 π x] Sin[π y]-4 π Cos[π y] Sin[π x] Sin[3 π x] Sin[3 π y]

VectorFieldPlot[{JX2[x,y],JY2[x,y]},{x,0,1},{y,0,1},PlotLabel->"Initial Probability Current: Example 2", LabelStyle->Directive[Bold, Red]]

VX2[x_,y_]:=((JX2[x,y])/(P2[x,y]+0.00000001))
VY2[x_,y_]:=((JY2[x,y])/(P2[x,y]+0.00000001))

VectorFieldPlot[{VX2[x,y],VY2[x,y]},{x,0,1},{y,0,1},PlotLabel->"Initial Velocity Field: Example 2", LabelStyle->Directive[Bold, Red]]

D[VY2[x,y],x]-D[VX2[x,y],y]//Simplify

Ω2[x_,y_]:=(7.106115168784338`*^-6 Cos[3 π x] Cos[3 π y] Sin[π x] Sin[π y]-7.895683520871487`*^-7 Cos[π x] Cos[π y] Sin[3 π x] Sin[3 π y])/(((1.`*^-8+2.` ((Sin[3 π x])^2) ((Sin[π y])^2)+2.` ((Sin[π x])^2) ((Sin[3 π y])^2)))^2)

Plot3D[Ω2[x,y],{x,0,1},{y,0,1},PlotPoints->100, PlotRange->{-100,100},Ticks->False,PlotLabel->"Initial Vorticity: Example 2", LabelStyle->Directive[Bold, Red]]

Plot3D[Arg[ψ2[x,y,0]],{x,0,1},{y,0,1},Ticks->{{},{},{-π,π}},PlotLabel->Style["Initial Phase: Example 2",Bold,Red]]

The singularities are positioned at {(1/3),(1/3)},  {(1/3),(2/3)}, {(2/3),(1/3)}, {(2/3),(2/3)}

Looking here only to the first of those, we have

Series[ψ2[(1/3)+dx,(1/3)+dy,t],{dx,0,1},{dy,0,1}]

so in first order

-3 i (Sqrt[(3/2)]) (e^(-10 i (π^2) t)) π dy-3 ((Sqrt[(3/2)]) (e^(-10 i (π^2) t)) π)dx

LocalPhase[dx_,dy_,t_]:=-3 I (Sqrt[(3/2)]) (E^(-10 I (π^2) t)) π dy-3 ((Sqrt[(3/2)]) (E^(-10 I (π^2) t)) π)dx

Manipulate[Plot3D[Arg[LocalPhase[dx,dy, (k/12)((2π)/(10(π^2)))]],{dx,((-1)/3),(2/3)},{dy,((-1)/3),(2/3)},Ticks->{{},{},{-π,π}},PlotLabel->Style["Time-dependence of phase near a singularity: Example 2",Bold,Red]],{k,0,11}]

THIRD EXAMPLE

MPF's final example uses coefficients borrowed from the 2nd example to assemble

ψ3[x_,y_,t_]:=(1/(Sqrt[2]))F31[x,y,t]+(I/(Sqrt[2]))F12[x,y,t]

ψ3[x,y,t]

ψ3star[x_,y_,t_]:=(Sqrt[2]) (E^(+10 I (π^2) t)) Sin[3 π x] Sin[π y]-I (Sqrt[2]) (E^(+5 I (π^2) t)) Sin[π x] Sin[2 π y]

ψ3star[x,y,t]ψ3[x,y,t]//ComplexExpand

P3[x_,y_,t_]:=2 ((Cos[10 (π^2) t])^2) ((Sin[3 π x])^2) ((Sin[π y])^2)+2 ((Sin[10 (π^2) t])^2) ((Sin[3 π x])^2) ((Sin[π y])^2)+4 Cos[10 (π^2) t] Sin[5 (π^2) t] Sin[π x] Sin[3 π x] Sin[π y] Sin[2 π y]-4 Cos[5 (π^2) t] Sin[10 (π^2) t] Sin[π x] Sin[3 π x] Sin[π y] Sin[2 π y]+2 ((Cos[5 (π^2) t])^2) ((Sin[π x])^2) ((Sin[2 π y])^2)+2 ((Sin[5 (π^2) t])^2) ((Sin[π x])^2) ((Sin[2 π y])^2)

Manipulate[Plot3D[P3[x,y,(k/12)((2π)/(5(π^2)))],{x,0,1},{y,0,1}, Ticks->False,PlotLabel->"Time-dependent Probability Density: Example 3", LabelStyle->Directive[Bold, Red]],{k,0,11}]

It is remarkable how briskly the distribution snaps from one fairly static shape to its mirror image. This snappy behavior is even more apparent when one looks to the probability current:

-I(ψ3star[x,y,t] D[ψ3[x,y,t],x]-ψ3[x,y,t] D[ψ3star[x,y,t],x])//ComplexExpand
-I(ψ3star[x,y,t] D[ψ3[x,y,t],y]-ψ3[x,y,t] D[ψ3star[x,y,t],y])//ComplexExpand

JX3[x_,y_,t_]:=-12 π Cos[5 (π^2) t] Cos[10 (π^2) t] Cos[3 π x] Sin[π x] Sin[π y] Sin[2 π y]-12 π Cos[3 π x] Sin[5 (π^2) t] Sin[10 (π^2) t] Sin[π x] Sin[π y] Sin[2 π y]+4 π Cos[5 (π^2) t] Cos[10 (π^2) t] Cos[π x] Sin[3 π x] Sin[π y] Sin[2 π y]+4 π Cos[π x] Sin[5 (π^2) t] Sin[10 (π^2) t] Sin[3 π x] Sin[π y] Sin[2 π y]

JY3[x_,y_,t_]:=
8 π Cos[5 (π^2) t] Cos[10 (π^2) t] Cos[2 π y] Sin[π x] Sin[3 π x] Sin[π y]+8 π Cos[2 π y] Sin[5 (π^2) t] Sin[10 (π^2) t] Sin[π x] Sin[3 π x] Sin[π y]-4 π Cos[5 (π^2) t] Cos[10 (π^2) t] Cos[π y] Sin[π x] Sin[3 π x] Sin[2 π y]-4 π Cos[π y] Sin[5 (π^2) t] Sin[10 (π^2) t] Sin[π x] Sin[3 π x] Sin[2 π y]

Manipulate[VectorFieldPlot[{JX3[x,y,(k/12)((2π)/(5(π^2)))],JY3[x,y,(k/12)((2π)/(5(π^2)))]},{x,0,1},{y,0,1},PlotLabel->"Time-dependent Probability Current: Example 3", LabelStyle->Directive[Bold, Red]],{k,0,11}]

I demonstrate that the vortices/singularities are actually located where they appear to be located:

P3[(1/3),(1/2),t]
P3[(2/3),(1/2),t]

VX3[x_,y_,t_]:=((JX3[x,y,t])/(P3[x,y,t]+0.00000001))
VY3[x_,y_,t_]:=((JY3[x,y,t])/(P3[x,y,t]+0.00000001))

Manipulate[VectorFieldPlot[{VX3[x,y,(k/12)((2π)/(5(π^2)))],VY3[x,y,(k/12)((2π)/(5(π^2)))]},{x,0,1},{y,0,1},PlotLabel->"Time-dependent Velocity Field: Example 3", LabelStyle->Directive[Bold, Red]],{k,0,11}]

D[VY3[x,y,0],x]-D[VX3[x,y,0],y]//Simplify

Ω3[x_,y_]:=(4.737410112522892`*^-6 Cos[3 π x] Cos[2 π y] Sin[π x] Sin[π y]-7.895683520871487`*^-7 Cos[π x] Cos[π y] Sin[3 π x] Sin[2 π y])/(((1.`*^-8+2.` ((Sin[3 π x])^2) ((Sin[π y])^2)+2.` ((Sin[π x])^2) ((Sin[2 π y])^2)))^2)

Plot3D[Ω3[x,y],{x,0,1},{y,0,1},PlotPoints->100, PlotRange->{-10,10}, Ticks->False]

Manipulate[Plot3D[Arg[ψ3[x,y,(k/12)((2π)/(5(π^2)))]],{x,0,1},{y,0,1},Ticks->{{},{},{-π,π}},PlotLabel->Style["Time-dependent Phase: Example 3",Bold,Red]],{k,0,11}]